package scu.maqiang.fes;

import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.transpose;
import static scu.maqiang.numeric.Constants.TetP3;
import static scu.maqiang.numeric.Constants.TetW3;


import java.util.Arrays;

import scu.maqiang.mesh.*;
import scu.maqiang.mesh.Mesh3T10;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES3T102 extends FES{

	public FES3T102(Mesh3T10 mesh) {
		super(mesh, 3);
		// TODO Auto-generated constructor stub
	}
	
	/**
	 * ��װ�ṹ����նȾ���, ���ϵ��Բ�����ϵ����������
	 * 
	 * @param coefFunc
	 *            ϵ��������2ά���飬������������ģ�������벴�ɱȺ���
	 * @param param
	 *            ϵ�������������
	 * @param tp
	 *            ��������
	 * @param A
	 *            ��װ�õ�������նȾ���
	 */
	public void assembleElasticStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES3T102::elementElasticityStiff, tp, A);
	}

	/**
	 * 
	 * @param constCoef
	 * @param tp
	 * @param A
	 */
	public void assembleElasticStiff(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES3T102::elementElasticityStiff, tp, A);
	}
	
	public void assembleElasticStiff(double[][][] matCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(matCoef, FES3T102::elementElasticityStiff, tp, A);
	}

	public void assembleElasticStiff(MatrixFunc coefFunc, double[][][] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES3T102::elementElasticityStiff, tp, A);
	}
	
	/**
	 * ���㵯�����ⵥԪ�նȾ���
	 * 
	 * @param xyz
	 *            3x2ά��Ԫ�������ÿ�б�ʾ����x��y����
	 * @param con
	 *            2x3ά���Գ���
	 * @param tp
	 * @param K
	 */
	public static void elementElasticityStiff(int[] ele, double[][] xyz, double[][] con, BVPType tp, double[][] K) {
		double[] N = new double[10];
		double[][] gradN = new double[10][3];
		double[] E = con[0];
		double[] Nu = con[1];
		double mE, mNu;
		double[][] B = null;
		double vol, weight;
		MVO.fill(K, 0.0);
		switch(tp) {
		case COMMON:
			for(int i = 0; i < 5; i++) {
				vol = CG.T10ShapeFunction(xyz, TetP3[i], N, gradN);
				mE = MVO.dot_product(E, N);
				mNu = MVO.dot_product(Nu, N);
				double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
				weight = vol * TetW3[i];
				B = new double[6][30];
				for(int ii = 0; ii < 10; ii++) {
					B[0][3 * ii]     = B[3][3 * ii + 1] = B[5][3 * ii + 2] = gradN[ii][0];
					B[1][3 * ii + 1] = B[3][3 * ii]     = B[4][3 * ii + 2] = gradN[ii][1];
					B[2][3 * ii + 2] = B[4][3 * ii + 1] = B[5][3 * ii]     = gradN[ii][2];
				}
				double[][] KK = matmul(transpose(B), matmul(D, B));
				for(int si = 0; si < 30; si++) {
					for(int sj = 0; sj < 30; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 5; i++) {
				vol = CG.T10ShapeFunction(xyz, TetP3[i], N, gradN);
				weight = vol * TetW3[i];
				B = new double[6][30];
				for(int ii = 0; ii < 10; ii++) {
					B[0][3 * ii]     = B[3][3 * ii + 1] = B[5][3 * ii + 2] = gradN[ii][0];
					B[1][3 * ii + 1] = B[3][3 * ii]     = B[4][3 * ii + 2] = gradN[ii][1];
					B[2][3 * ii + 2] = B[4][3 * ii + 1] = B[5][3 * ii]     = gradN[ii][2];
				}
				double[][] KK = matmul(transpose(B), matmul(con, B));
				for(int si = 0; si < 30; si++) {
					for(int sj = 0; sj < 30; sj++) {
						K[si][sj] += weight * KK[si][sj];
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleElasticMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES3T102::elementElasticityMass, tp, A);
	}

	public void assembleElasticMass(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES3T102::elementElasticityMass, tp, A);
	}

	public static void elementElasticityMass(int[] ele, double[][] xy, double[][] rho, BVPType tp, double[][] K) {
//		double[] p = rho[0];
//		double[][] grad = new double[3][2];
//		double area = CG.gradient2DT3(xy, grad);
//		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
//		if (tp == BVPType.CONSTITUTE_MATRIX_AXIS || tp == BVPType.AXIS_LUMP_MATRIX) {
//			MVO.mul(p, x, p);
//		}
//		double N11 = (3 * p[0] + p[1] + p[2]) / 30 * area;
//		double N22 = (p[0] + 3 * p[1] + p[2]) / 30 * area;
//		double N33 = (p[0] + p[1] + 3 * p[2]) / 30 * area;
//		double N12 = (2 * p[0] + 2 * p[1] + p[2]) / 60 * area;
//		double N13 = (2 * p[0] + p[1] + 2 * p[2]) / 60 * area;
//		double N23 = (p[0] + 2 * p[1] + 2 * p[2]) / 60 * area;
//		MVO.fill(K, 0.0);
//		switch (tp) {
//		case CONSIST_MATRIX:
//		case AXIS_CONSIST_MATRIX:
//			K[0][0] = K[1][1] = N11;
//			K[0][2] = K[1][3] = N12;
//			K[0][4] = K[1][5] = N13;
//
//			K[2][0] = K[3][1] = N12;
//			K[2][2] = K[3][3] = N22;
//			K[2][4] = K[3][5] = N23;
//
//			K[4][0] = K[5][1] = N13;
//			K[4][2] = K[5][3] = N23;
//			K[4][4] = K[5][5] = N33;
//
//			break;
//		case LUMP_MATRIX:
//		case AXIS_LUMP_MATRIX:
//			K[0][0] = K[1][1] = (2 * p[0] + p[1] + p[2]) / 12 * area;
//			K[2][2] = K[3][3] = (p[0] + 2 * p[1] + p[2]) / 12 * area;
//			K[4][4] = K[5][5] = (p[0] + p[1] + 2 * p[2]) / 12 * area;
//			break;
//		default:
//			break;
//		}
		// System.out.println(MVO.toString(K));
		// System.exit(0);
	}

	public void assembleVolumeForce(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES3T102::elementVolumeForce, tp, RHS);
	}

	public void assembleVolumeForce(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES3T102::elementVolumeForce, tp, RHS);
	}

	public static void elementVolumeForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] N = new double[10];
		double[][] gradN = new double[10][3];
		double mfu, mfv, mfw;
		double vol, weight;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 5; i++) {
				vol = CG.T10ShapeFunction(xy, TetP3[i], N, gradN);
				mfu = MVO.dot_product(f[0], N);
				mfv = MVO.dot_product(f[1], N);
				mfw = MVO.dot_product(f[2], N);
				weight = vol  * TetW3[i];
				for (int si = 0; si < 10; si++) {
					eleF[3 * si]     += N[si] * weight * mfu;
					eleF[3 * si + 1] += N[si] * weight * mfv;
					eleF[3 * si + 2] += N[si] * weight * mfw;
				}				
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleSurfaceForce(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES3T102::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public void assembleSurfaceForce(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES3T102::elemenSurfaceForce, tp, RHS, bdLabel);
	}
	
	public void assembleSurfaceForce(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES3T102::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public static void elemenSurfaceForce(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double area = CG.area3DT3(xyz);
//		double[][] B = new double[6][6];
//		double[] fuv = {f[0][0], f[1][0], f[0][1], f[1][1], f[0][2], f[1][2]};
//		B[0][0] = B[1][1] = B[2][2] = B[3][3] = 2.0 * length / 15;
//		B[0][2] = B[2][0] = B[1][3] = B[3][1] = -length / 30.0;
//		B[0][4] = B[1][5] = B[2][4] = B[3][5] = length / 15.0;
//		B[4][0] = B[4][2] = B[5][1] = B[5][3] = length / 15.0;
//		B[4][4] = B[5][5] = 8.0 * length / 15.0;
//		MVO.matmul(B, fuv, ele);
		
		double[] gu = f[0];
		double[] gv = f[1];
		double[] gw = f[2];
		double[][] L = {{area/30,   -area/180, -area/180, 0, -area/45, 0}, 
				        {-area/180, area/30,   -area/180, 0,        0, -area/45},
				        {-area/180, -area/180, area/30, -area/45, 0, 0},
				        {0, 0, -area/45, 8*area/45, 4*area/45, 4*area/45},
				        {-area/45, 0, 0, 4*area/45, 8*area/45, 4*area/45},
				        {0, -area/45, 0, 4*area/45, 4*area/45, 8*area/45}};
		switch (tp) {
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 6; i++) {
				eleF[3 * i] = MVO.dot_product(L[i], gu);
				eleF[3 * i + 1] = MVO.dot_product(L[i], gv);
				eleF[3 * i + 2] = MVO.dot_product(L[i], gw);
			}
//			ele[0] = (4 * gu[0] - gu[1] + 2 * gu[2]) * area / 30;
//			ele[1] = (4 * gv[0] - gv[1] + 2 * gv[2]) * area / 30;
//			ele[2] = (-gu[0] + 4 * gu[1] + 2 * gu[2]) * area / 30;
//			ele[3] = (-gv[0] + 4 * gv[1] + 2 * gv[2]) * area / 30;
//			ele[4] = (gu[0] + gu[1] + 8 * gu[2]) * area / 15;
//			ele[5] = (gv[0] + gv[1] + 8 * gv[2]) * area / 15;
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		
	}

	public void assemblePressure(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES3T102::elementPressure, tp, RHS, bdLabel);
	}

	public void assemblePressure(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES3T102::elementPressure, tp, RHS, bdLabel);
	}
	
	public void assemblePressure(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES3T102::elementPressure, tp, RHS, bdLabel);
	}
	
	public static void elementPressure(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] pressure = f[0];
		double[] n = CG.unitOutNormal3DT3(xyz);
//		f[0][0] = -pressure[0] * n[0];
//		f[0][1] = -pressure[1] * n[0];
//		f[0][2] = -pressure[2] * n[0];
//		f[0][3] = -pressure[3] * n[0];
//		f[0][4] = -pressure[4] * n[0];
//		f[0][5] = -pressure[5] * n[0];
//		f[1][0] = -pressure[0] * n[1];
//		f[1][1] = -pressure[1] * n[1];
//		f[1][2] = -pressure[2] * n[1];
//		f[1][3] = -pressure[3] * n[1];
//		f[1][4] = -pressure[4] * n[1];
//		f[1][5] = -pressure[5] * n[1];
//		f[2][0] = -pressure[0] * n[2];
//		f[2][1] = -pressure[1] * n[2];
//		f[2][2] = -pressure[2] * n[2];
//		f[2][3] = -pressure[3] * n[2];
//		f[2][4] = -pressure[4] * n[2];
//		f[2][5] = -pressure[5] * n[2];
		for(int i = 0; i < 3; i++) {
			for(int j = 0; j < 6; j++) {
				f[i][j] = -pressure[j] * n[i];
			}
		}
		elemenSurfaceForce(ele, xyz, f, tp, eleF);
	}	

	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES3T102::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES3T102::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
//		double[] beta = coef[0];
//		double[] temperature = coef[1];
//		double[] x = MVO.col(coord, 0);
//		double[][] grad = new double[3][2];
//		double area = CG.gradient2DT3(coord, grad);
//		Arrays.fill(eleLoad, 0.0);
//		switch (tp) {
//		case PLANESTRAIN:
//		case PLANESTRESS:
//		case CONSTITUTE_MATRIX_COMMON:
//			double c = (MVO.sum(beta) * MVO.sum(temperature) + MVO.dot_product(beta, temperature)) * area / 12.0;
//			for (int i = 0; i < 3; i++) {
//				eleLoad[2 * i] = c * grad[i][0];
//				eleLoad[2 * i + 1] = c * grad[i][1];
//			}
//			break;
//		case AXISYMMETRIC:
//		case CONSTITUTE_MATRIX_AXIS:
//			for (int i = 0; i < 4; i++) {
//				double[] L = Constants.TP3[i];
//				double wxarea = MVO.dot_product(beta, L) * MVO.dot_product(temperature, L) * TW3[i];
//				for (int j = 0; j < 3; j++) {
//					eleLoad[2 * j] += (grad[j][0] * MVO.dot_product(x, L) + L[j]) * wxarea;
//					eleLoad[2 * j + 1] += grad[j][1] * MVO.dot_product(x, L) * area;
//				}
//			}
//			break;
//		default:
//			throw new IllegalArgumentException("Wrong BVPType!");
//		}
		// System.out.println(Arrays.toString(eleLoad));

	}

	public void computeNodalStrain(double[] sol, double[][] strain) {
//		int len = sol.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		int n = len / 2;
//		ParamCheck.checkEqual(n, strain[0].length);
//		ParamCheck.checkEqual(n, strain[1].length);
//		ParamCheck.checkEqual(n, strain[2].length);
//		double[] angle = new double[3];
//		double[][] coord = new double[3][2];
//		double[][] grad = new double[3][2];
//		double[] angleCount = new double[this.GetNdof()];
//		double dxU, dyU, dxV, dyV;
//		Element ele = null;
//		double[] eleU = new double[3];
//		double[] eleV = new double[3];
//		for (int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.gradient2DT3(coord, grad);
//			CG.angle2DT3(coord, angle);
//			for (int j = 0; j < 3; j++) {
//				eleU[j] = sol[2 * ele.idx[j]];
//				eleV[j] = sol[2 * ele.idx[j] + 1];
//			}
//			double[] dx = MVO.col(grad, 0);
//			double[] dy = MVO.col(grad, 1);
//			dxU = MVO.dot_product(eleU, dx);
//			dyU = MVO.dot_product(eleU, dy);
//			dxV = MVO.dot_product(eleV, dx);
//			dyV = MVO.dot_product(eleV, dy);
//			for (int j = 0; j < 3; j++) {
//				strain[0][ele.idx[j]] += dxU * angle[j];
//				strain[1][ele.idx[j]] += dyV * angle[j];
//				strain[2][ele.idx[j]] += (dxV + dyU) * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for (int i = 0; i < n; i++) {
//			strain[0][i] /= angleCount[i];
//			strain[1][i] /= angleCount[i];
//			strain[2][i] /= angleCount[i];
//		}
		
//		int len = u.length;
//		ParamCheck.checkEqual(len, this.ndof);
//
//		double[] angle = new double[6];
//		double[][] coord = new double[6][2];
//		double[][] grad = new double[3][2];
//		double[][] eleGradU = new double[6][2];
//		double[] eleSDU = new double[4];
//		double[] angleCount = new double[this.GetNdof()];
//		Element ele = null;
//		double[] eleU = new double[6];
//		for(int i = 0, n = mesh.getNt(); i < n; i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = u[ele.idx[j]];
//			}
//			CG.gradient2DT3(coord, grad);
//			CG.gradientValue2DT6(grad, eleU, eleGradU);
//			CG.secondDerivative2DT6(grad, eleU, eleSDU);
//			CG.angle2DT6(coord, angle);
//			for(int j = 0; j < 6; j++) {
//				fdU[0][ele.idx[j]] += eleGradU[j][0] * angle[j];
//				fdU[1][ele.idx[j]] += eleGradU[j][1] * angle[j];
//				sdU[0][ele.idx[j]] += eleSDU[0] * angle[j];
//				sdU[1][ele.idx[j]] += eleSDU[1] * angle[j];
//				sdU[2][ele.idx[j]] += eleSDU[2] * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for(int i = 0, n = mesh.getNv(); i < n; i++) {
//			fdU[0][i] /= angleCount[i];
//			fdU[1][i] /= angleCount[i];
//			sdU[0][i] /= angleCount[i];
//			sdU[1][i] /= angleCount[i];
//			sdU[2][i] /= angleCount[i];
//		}
//		
//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[0][ele.idx[j]];
//			}
//			sdU[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//		}
	}


	/**
	 * 
	 * @param uvw
	 * @param fdUVW
	 * @param sdUVW
	 */
	public void computeFirstAndSecondDerivative(double[][] uvw, double[][] fdUVW, double[][] sdUVW) {
		int len = uvw[0].length;
		ParamCheck.checkEqual(len * 3, this.ndof);
		double[][] coord = new double[10][3];
		double[][] grad = new double[4][3];
		double[][] eleGradU = new double[10][3];
		double[][] eleGradV = new double[10][3];
		double[][] eleGradW = new double[10][3];
		double[] eleSDU = new double[10];
		double[] eleSDV = new double[10];
		double[] eleSDW = new double[10];
		double[] volumeCount = new double[this.GetNdof()];
		double[] eleU = new double[10];
		double[] eleV = new double[10];
		double[] eleW = new double[10];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 10; j++) {
				eleU[j] = uvw[0][ele[j]];
				eleV[j] = uvw[1][ele[j]];
				eleW[j] = uvw[2][ele[j]];
			}
			double volume = CG.gradient3DT4(coord, grad);
			CG.gradientValue3DT10(grad, eleU, eleGradU);
			CG.gradientValue3DT10(grad, eleV, eleGradV);
			CG.gradientValue3DT10(grad, eleW, eleGradW);
			CG.secondDerivative3DT10(grad, eleU, eleSDU);
			CG.secondDerivative3DT10(grad, eleV, eleSDV);
			CG.secondDerivative3DT10(grad, eleW, eleSDW);
			for(int j = 0; j < 10; j++) {
				fdUVW[0][ele[j]] += eleGradU[j][0] * volume;
				fdUVW[1][ele[j]] += eleGradU[j][1] * volume;
				fdUVW[2][ele[j]] += eleGradU[j][2] * volume;
				
				fdUVW[3][ele[j]] += eleGradV[j][0] * volume;
				fdUVW[4][ele[j]] += eleGradV[j][1] * volume;
				fdUVW[5][ele[j]] += eleGradV[j][2] * volume;
				
				fdUVW[6][ele[j]] += eleGradW[j][0] * volume;
				fdUVW[7][ele[j]] += eleGradW[j][1] * volume;
				fdUVW[8][ele[j]] += eleGradW[j][2] * volume;
				
				sdUVW[0][ele[j]] += eleSDU[0] * volume;
				sdUVW[1][ele[j]] += eleSDU[1] * volume;
				sdUVW[2][ele[j]] += eleSDU[2] * volume;
				sdUVW[3][ele[j]] += eleSDU[3] * volume;
				sdUVW[4][ele[j]] += eleSDU[4] * volume;
				sdUVW[5][ele[j]] += eleSDU[5] * volume;
				
				sdUVW[6][ele[j]] += eleSDV[0] * volume;
				sdUVW[7][ele[j]] += eleSDV[1] * volume;
				sdUVW[8][ele[j]] += eleSDV[2] * volume;
				sdUVW[9][ele[j]] += eleSDV[3] * volume;
				sdUVW[10][ele[j]] += eleSDV[4] * volume;
				sdUVW[11][ele[j]] += eleSDV[5] * volume;
				
				sdUVW[12][ele[j]] += eleSDW[0] * volume;
				sdUVW[13][ele[j]] += eleSDW[1] * volume;
				sdUVW[14][ele[j]] += eleSDW[2] * volume;
				sdUVW[15][ele[j]] += eleSDW[3] * volume;
				sdUVW[16][ele[j]] += eleSDW[4] * volume;
				sdUVW[17][ele[j]] += eleSDW[5] * volume;
				volumeCount[ele[j]] += volume;
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			fdUVW[0][i] /= volumeCount[i];
			fdUVW[1][i] /= volumeCount[i];
			fdUVW[2][i] /= volumeCount[i];
			fdUVW[3][i] /= volumeCount[i];
			fdUVW[4][i] /= volumeCount[i];
			fdUVW[5][i] /= volumeCount[i];
			fdUVW[6][i] /= volumeCount[i];
			fdUVW[7][i] /= volumeCount[i];
			fdUVW[8][i] /= volumeCount[i];
			sdUVW[0][i] /= volumeCount[i];
			sdUVW[1][i] /= volumeCount[i];
			sdUVW[2][i] /= volumeCount[i];
			sdUVW[3][i] /= volumeCount[i];
			sdUVW[4][i] /= volumeCount[i];
			sdUVW[5][i] /= volumeCount[i];
			sdUVW[6][i] /= volumeCount[i];
			sdUVW[7][i] /= volumeCount[i];
			sdUVW[8][i] /= volumeCount[i];
			sdUVW[9][i] /= volumeCount[i];
			sdUVW[10][i] /= volumeCount[i];
			sdUVW[11][i] /= volumeCount[i];
			sdUVW[12][i] /= volumeCount[i];
			sdUVW[13][i] /= volumeCount[i];
			sdUVW[14][i] /= volumeCount[i];
			sdUVW[15][i] /= volumeCount[i];
			sdUVW[16][i] /= volumeCount[i];
			sdUVW[17][i] /= volumeCount[i];
		}
		
		for(int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			for(int k = 0; k < 18; k++) {
				for(int j = 0; j < 10; j++) {
					eleU[j] = sdUVW[k][ele[j]];
				}
				sdUVW[k][ele[4]] = 0.5 * (eleU[0] + eleU[1]);
				sdUVW[k][ele[5]] = 0.5 * (eleU[1] + eleU[2]);
				sdUVW[k][ele[6]] = 0.5 * (eleU[2] + eleU[3]);
				sdUVW[k][ele[7]] = 0.5 * (eleU[0] + eleU[3]);
				sdUVW[k][ele[8]] = 0.5 * (eleU[1] + eleU[3]);
				sdUVW[k][ele[9]] = 0.5 * (eleU[2] + eleU[3]);			
			}
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[6]] = 0.5 * (eleU[2] + eleU[3]);
//			sdU[1][ele.idx[7]] = 0.5 * (eleU[0] + eleU[3]);
//			sdU[1][ele.idx[8]] = 0.5 * (eleU[1] + eleU[3]);
//			sdU[1][ele.idx[9]] = 0.5 * (eleU[2] + eleU[3]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[6]] = 0.5 * (eleU[2] + eleU[3]);
//			sdU[2][ele.idx[7]] = 0.5 * (eleU[0] + eleU[3]);
//			sdU[2][ele.idx[8]] = 0.5 * (eleU[1] + eleU[3]);
//			sdU[2][ele.idx[9]] = 0.5 * (eleU[2] + eleU[3]);
		}
			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[1][ele.idx[j]];
//			}
//			sdUV[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[2][ele.idx[j]];
//			}
//			sdUV[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[0][ele.idx[j]];
//			}
//			sdUV[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[1][ele.idx[j]];
//			}
//			sdUV[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdUV[2][ele.idx[j]];
//			}
//			sdUV[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdUV[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdUV[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
		
	}

}
